

clear
set more off
cap log close

***************************************************************************************************
* Program: figureS3.do
* Purpose: Produce figure S3
***************************************************************************************************

************************************************************************
** set macros (directory)
************************************************************************
global prison_dir "~/Dropbox/Prison_Covid/PNAS_Nexus_Replication"
global raw_data "${prison_dir}/raw_data"
global intermediate_data "${prison_dir}/intermediate_data"
global code "${prison_dir}/code"
global analysis_data "${prison_dir}/analysis_data"
global output "${prison_dir}/output"

* Variables 
global demo "pct_nh_blk_alone_acs_15_19 pct_nh_white_alone_acs_15_19 pct_hisp_latino_acs_15_19 pct_college_acs_15_19 med_hhc_inc_acs_15_19 pct_pop_65plus_acs_15_19 log_tot_pop_acs_15_19"
global cubic "ZCTA_max_cases_p100k_rmpri5 ZCTA_max_cases_p100k_rmpri5_sq ZCTA_max_cases_p100k_rmpri5_cu"

global covidXwhite "cum_cases5_X_pct_white cases_5_X_pct_white cases_sq_5_X_pct_white cases_cu_5_X_pct_white"
global covidXinc "cum_cases5_X_inc cases_5_X_inc cases_sq_5_X_inc cases_cu_5_X_inc"	


***************************************
* Main Code
***************************************
use "${analysis_data}/San_Quentin_main_analysis_data.dta", clear 
	
***************************************
* get pscore and wgt under different matching schemes * 
***************************************

drdid ZCTA_max_cases_p100k_rmpri $cubic cum_max_cases_p100k_rmpri5 log_tot_pop_acs_15_19 if inlist(month, 5, 6), ///
	ivar(ZCTA_contact) time(month) tr(treated) dripw
drdid_predict ps_covid, pscore
drdid_predict wgt_covid_temp 

drdid ZCTA_max_cases_p100k_rmpri $cubic cum_max_cases_p100k_rmpri5 $demo if inlist(month, 5, 6), ///
	ivar(ZCTA_contact) time(month) tr(treated) dripw
drdid_predict ps_covid_demo, pscore
drdid_predict wgt_covid_demo_temp


gegen wgt_covid = max(wgt_covid_temp), by(ZCTA_contact)
gegen wgt_covid_demo = max(wgt_covid_demo_temp), by(ZCTA_contact)

drop *_temp

***************************************
* density plot of pscore 
***************************************
local name = "San Quentin"
tw (kdensity ps_covid if treated == 1 & month == 5, lcolor(red*.5)  fcolor(none)) ///
   (kdensity ps_covid if treated == 0 & month == 5, ///
		 lcolor(blue*.5)  fcolor(none) ///
		legend(order(1 "Zip Codes Connected to `name' State Prison" ///
		2 "Zip Codes Unconnected to `name' State Prison") row(2) position(6)) ///
		xtitle("Propensity Score (Unweighted)") ytitle("Density") ///
		title("Panel A: Matched on COVID-19 Variables", size(medsmall))) 
graph save "${output}/kdensity_csdid_pscore_covid_SQ.gph", replace 
		
tw (kdensity ps_covid if treated == 1 & month == 5, lcolor(red*.5)  fcolor(none)) ///
   (kdensity ps_covid if treated == 0 & month == 5 [w = wgt_covid], ///
		 lcolor(blue*.5)  fcolor(none) ///
		legend(order(1 "Zip Codes Connected to `name' State Prison" ///
		2 "Zip Codes Unconnected to `name' State Prison") row(2) position(6)) ///
		xtitle("Propensity Score (Weighted)") ytitle("Density") ///
		title("Panel B: Matched on COVID-19 Variables", size(medsmall))) 
graph save "${output}/kdensity_csdid_pscore_covid_wgt_SQ.gph", replace 

tw (kdensity ps_covid_demo if treated == 1 & month == 5, lcolor(red*.5)  fcolor(none)) ///
   (kdensity ps_covid_demo if treated == 0 & month == 5, ///
		 lcolor(blue*.5)  fcolor(none) ///
		legend(order(1 "Zip Codes Connected to `name' State Prison" ///
		2 "Zip Codes Unconnected to `name' State Prison") row(2) position(6)) ///
		xtitle("Propensity Score (Unweighted)") ytitle("Density") ///
		title("Panel C: Matched on COVID-19 Variables + Demographics", size(medsmall))) 
graph save "${output}/kdensity_csdid_pscore_covid_demo_SQ.gph", replace 
		
tw (kdensity ps_covid_demo if treated == 1 & month == 5, lcolor(red*.5)  fcolor(none)) ///
   (kdensity ps_covid_demo if treated == 0 & month == 5 [w = wgt_covid_demo], ///
		 lcolor(blue*.5)  fcolor(none) ///
		legend(order(1 "Zip Codes Connected to `name' State Prison" ///
		2 "Zip Codes Unconnected to `name' State Prison") row(2) position(6)) ///
		xtitle("Propensity Score (Weighted)") ytitle("Density") ///
		title("Panel D: Matched on COVID-19 Variables + Demographics", size(medsmall))) 
graph save "${output}/kdensity_csdid_pscore_covid_demo_wgt_SQ.gph", replace 
		
** combine graphs:
cd "${output}"
graph combine ///
	kdensity_csdid_pscore_covid_SQ.gph kdensity_csdid_pscore_covid_wgt_SQ.gph ///
	kdensity_csdid_pscore_covid_demo_SQ.gph kdensity_csdid_pscore_covid_demo_wgt_SQ.gph, ///
	xcommon  
graph export "${output}/figures3.png", replace

rm kdensity_csdid_pscore_covid_SQ.gph
rm kdensity_csdid_pscore_covid_wgt_SQ.gph
rm kdensity_csdid_pscore_covid_demo_SQ.gph
rm kdensity_csdid_pscore_covid_demo_wgt_SQ.gph

